#########################R-script for bookchapter: Individual learning behaviour in collaborative processes of innovation ###################################
######################### Written by Vidar Stevens - vidar.stevens@uantwerpen.be #################################################

##################Strategy of analysis##########################
#To perform the ERGM analyses I followed the analysis and model selection strategy that has been proposed by Goodreau (2007)
#For each dimension of the learning variable (i.e. construction, co-construction, and constructive conflict) four, and sometimes five models, will be tested.
#A first model is always the so-called Bernoulli/Erd�s-R�nyi model, which with only a single term (i.e. edges) captures the density of the network under study (Goodreau, 2007:239).
#Second, I include the selected endogenous terms in the analysis. Within ERGM studies this is considered as a 'reduced homogeneous realization-dependent model' (idem: 240).
#Third, I add the exogenous nodal covariates (i.e. amount of years representative, top-priority Minister, type of organization, personal meaningfulness, and chairman) to the model with the endogenous terms. 
#Fourth, I analyse a model that includes the endogenous terms as well as the exogenous dyadic covariate effects (i.e. belief homophily, trust, and actor importance).
#If necessary, I test a fifth model - a more optimized model of the dyadic covariate effects model - to ensure a better fit of the (final) explanatory model. 
#Hence, for each dimension of the dependent variable I search for four, or sometimes five, 'best fitting models' - and from this sample of models I select the model that is most capable in explaining the learning interaction patters in the empirical data. 

##################Model-selection and goodness of fit criteria######################
#Inspired by the work of Goodreau (2007:238), I use three general approaches to examine the goodness of fit of the models. 
#First, I check for degeneracy and model convergence. A minimum requirement for a model to fit well is for estimation of parameters to converge on finite parameter values. It must also be non-degenerate, that is, not place all of its probability mass on a few networks entirely unlike the observed network, such as a full or empty network.
#Second, I compare the Akaike Information Criterion (AIC) between models. Models that exhibit dyadic independence can be fit with standard logistic regressions, which yield a likelihood measure for the model. Those models that are dyad dependent must be fit with MCMC, which also yields an estimate of the likelihood. I use the given likelihood to calculate AIC in order to compare models; with a lower AIC implying a significant increase in model fit. 
#Third, I view the goodness of fit by plotting the model against higher order statistics. This approach is described in detail by Hunter et al. (2008). The logic entails generating new networks according to the probability distribution implied by the fit model. Because the normalizing constant is still present in the fit models, this must be done using the same MCMC approach employed during the estimation procedure. A statistic of interest is then calculated on the original network and on the set of networks generated from the model, and these are plotted for comparison. If the original network is inconsistent with the networks generated from the model, this suggests that the structure of the network differs from those predicted by the model, and the model is not well fit. Multiple statistics can be compared visually to provide detailed information about the systematic ways in which the data and the model predictions differ.
#For this third approach, the higher order network statistics I specifically use to compare the selected models to include: in-degree, out-degree, edge-wise shared partners, dyad-wise shared partners, and the minimum geodesic distance. Each of these higher-order network statistics are plotted on the log-odds scale for the sake of greater visibility across the covered range of values. 

###################Neccesary packages to use the ERGM method########################
#Install the neccesary packages, load the library, set work directory and read in the csv.file
install.packages("statnet","ergm","sna","boot")
install.packages("latticeExtra")
library("statnet")
library("sna")
library("boot")
library("ergm")
library("latticeExtra")

###############R-codes per dimension#####################
#Our dependent variable of learning has 3 seperate dimensions: construction, co-construction and constructive conflict.
#These are derived from different scales
#For the analysis of each of these dimensions we will first do the preparations.
#We start with the dimension of construction.

###############construction (preperation)################
#setwd, load table, turn it into adjacency-matrix.
#In the survey we asked each respondent to indicate up to 5 representatives that showed signs of 'construction'learning behaviour to them.
#Hence we transpose the excel matrix, to get a table representating out-going ties.
setwd("C:/Users/Vstevens/Box Sync/PhD/Vlaamse contacten/Ambtelijk Forum BRV/Survey-responses")
construction=read.csv(("Construction_transponeren.csv"), header=T, row.names=1, sep=';')

#turn the data-files into a matrix-object
m=as.matrix(construction)

#turn the matrix-object into a network adjacency matrix and plot the graph
net=network(m, matrix.type="adjacency", directed=TRUE)
plot(net,displaylabels=T)
summary(net)
network.dyadcount(net)
network.edgecount(net)
network.size(net)
network.density(net)
plot(net,displaylabels=T, mode="circle")

#add datasetAF  
datasetAF<- read.table("datasetAF.csv", header=T, row.names=1, sep=";", stringsAsFactors=FALSE)

#convert attributes of datasetAF to vectors
years_experience<-datasetAF[,1]
years_member_AF <- datasetAF[,2]
toppriority_minister <- datasetAF[,17]
toppriority_administrativeleader <- datasetAF[,18]
controll_minister<-datasetAF[,19]
controll_administrativeleader <-datasetAF[,20]
experienced_autonomy <-datasetAF[,21]
societalvalue <- datasetAF[,22]
personalvalue <- datasetAF[,23]
chairman1<- datasetAF[,24]
projectleader1 <- datasetAF[,25]
agency1 <- datasetAF[,26]
ruimteclaim1<- datasetAF[,27]
hardeb<- datasetAF[,28]
zachteb<- datasetAF[,29]
kwaltb<- datasetAF[,30]
logistieksysteem1<- datasetAF[,31]
representative1<- datasetAF[,32]

#attach vectors to the network net
set.vertex.attribute(net,"yearsexperience",years_experience)
set.vertex.attribute(net,"yearsrepresentative",years_member_AF)
set.vertex.attribute(net,"controllminister",controll_minister)
set.vertex.attribute(net,"controlladministrativeleader",controll_administrativeleader)
set.vertex.attribute(net,"persvalue",personalvalue)
set.vertex.attribute(net,"societvalue",societalvalue)
set.vertex.attribute(net,"toppminister",toppriority_minister)
set.vertex.attribute(net,"toppadml",toppriority_administrativeleader)
set.vertex.attribute(net,"autonomy",experienced_autonomy)
set.vertex.attribute(net,"chairman",chairman1) #with this vertex attribute, we identify the chairman in the network
set.vertex.attribute(net,"agency",agency1) 
set.vertex.attribute(net,"ruimteclaim",ruimteclaim1) 
set.vertex.attribute(net,"hardebelangen",hardeb) 
set.vertex.attribute(net,"zachtebelangen",zachteb) 
set.vertex.attribute(net,"kwaltbelangen",kwaltb) 
set.vertex.attribute(net,"logistieksysteem",logistieksysteem1) 
set.vertex.attribute(net,"projectleader",projectleader1) 
set.vertex.attribute(net,"representative",representative1) 

#add covariates beliefhomophily, trust and earlier billateral contact patterns
beliefhomophily<- as.matrix(read.table("beliefhomo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
trust<- as.matrix(read.table("trust.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
earlierbillateralcontact<- as.matrix(read.table("earliercontact.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
actorimportance<- as.matrix(read.table("actor_impo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))

diag(beliefhomophily)<- 0
diag(trust)<- 0
diag(earlierbillateralcontact)<- 0
diag(actorimportance)<- 0

###############end of preperation dimension of construction#####################

###############co-construction(preparation)###########################
#In the survey we asked each respondent to indicate up to 5 representatives that showed signs of 'co-construction'learning behaviour to them.
#Hence we transpose the excel matrix, to get a table representating out-going ties.
setwd("C:/Users/Vstevens/Box Sync/PhD/Vlaamse contacten/Ambtelijk Forum BRV/Survey-responses")
coconstruction=read.csv(("coconstruction_transponeren.csv"), header=T, row.names=1, sep=';')

#turn the data-files into a matrix-object
mm=as.matrix(coconstruction)

#turn the matrix-object into a network adjacency matrix and plot the graph
netmm=network(mm, matrix.type="adjacency", directed=TRUE)
plot(netmm,displaylabels=TRUE)
summary(netmm)
network.dyadcount(netmm)
network.edgecount(netmm)
network.size(netmm)
network.density(netmm)
plot(netmm,displaylabels=T, mode="circle")

#add datasetAF  
datasetAF<- read.table("datasetAF.csv", header=T, row.names=1, sep=";", stringsAsFactors=FALSE)

#convert attributes of datasetAF to vectors
years_experience<-datasetAF[,1]
years_member_AF <- datasetAF[,2]
toppriority_minister <- datasetAF[,17]
toppriority_administrativeleader <- datasetAF[,18]
controll_minister<-datasetAF[,19]
controll_administrativeleader <-datasetAF[,20]
experienced_autonomy <-datasetAF[,21]
societalvalue <- datasetAF[,22]
personalvalue <- datasetAF[,23]
chairman1<- datasetAF[,24]
projectleader1 <- datasetAF[,25]
agency1 <- datasetAF[,26]
ruimteclaim1<- datasetAF[,27]
hardeb<- datasetAF[,28]
zachteb<- datasetAF[,29]
kwaltb<- datasetAF[,30]
logistieksysteem1<- datasetAF[,31]

#attach vectors to the network netmm
set.vertex.attribute(netmm,"yearsexperience",years_experience)
set.vertex.attribute(netmm,"yearsrepresentative",years_member_AF)
set.vertex.attribute(netmm,"controllminister",controll_minister)
set.vertex.attribute(netmm,"controlladministrativeleader",controll_administrativeleader)
set.vertex.attribute(netmm,"persvalue",personalvalue)
set.vertex.attribute(netmm,"societvalue",societalvalue)
set.vertex.attribute(netmm,"toppminister",toppriority_minister)
set.vertex.attribute(netmm,"toppadml",toppriority_administrativeleader)
set.vertex.attribute(netmm,"autonomy",experienced_autonomy)
set.vertex.attribute(netmm,"chairman",chairman1) #with this vertex attribute, we identify the chairman in the network
set.vertex.attribute(netmm,"projectleader",projectleader1)
set.vertex.attribute(netmm,"agency",agency1) 
set.vertex.attribute(netmm,"ruimteclaim",ruimteclaim1) 
set.vertex.attribute(netmm,"hardebelangen",hardeb) 
set.vertex.attribute(netmm,"zachtebelangen",zachteb) 
set.vertex.attribute(netmm,"kwaltbelangen",kwaltb) 
set.vertex.attribute(netmm,"logistieksysteem",logistieksysteem1) 



#add covariates beliefhomophily, trust and earlier billateral contact patterns
beliefhomophily<- as.matrix(read.table("beliefhomo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
trust<- as.matrix(read.table("trust.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
earlierbillateralcontact<- as.matrix(read.table("earliercontact.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
actorimportance<- as.matrix(read.table("actor_impo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))

diag(beliefhomophily)<- 0
diag(trust)<- 0
diag(earlierbillateralcontact)<- 0
diag(actorimportance)<- 0

#######################end of preperation dimension of co-construction#################################

############################constructive conflict (preparation)##################################
#In the survey we asked each respondent to indicate up to 5 representatives that showed signs of 'constructive conflict'learning behaviour to them.
#Hence we transpose the excel matrix, to get a table representating out-going ties.
setwd("C:/Users/Vstevens/Box Sync/PhD/Vlaamse contacten/Ambtelijk Forum BRV/Survey-responses")
constructiveconflict=read.csv(("constructiveconflict_transponeren.csv"), header=T, row.names=1, sep=';')

#turn the data-files into a matrix-object
mmm=as.matrix(constructiveconflict)

#turn the matrix-object into a network adjacency matrix and plot the graph
netmmm=network(mmm, matrix.type="adjacency", directed=TRUE)
plot(netmmm,displaylabels=TRUE)
summary(netmmm)
network.dyadcount(netmmm)
network.edgecount(netmmm)
network.size(netmmm)
network.density(netmmm)
plot(netmmm,displaylabels=F, mode="circle")

#add datasetAF  
datasetAF<- read.table("datasetAF.csv", header=T, row.names=1, sep=";", stringsAsFactors=FALSE)

#convert attributes of datasetAF to vectors#
years_experience<-datasetAF[,1]
years_member_AF <- datasetAF[,2]
toppriority_minister <- datasetAF[,17]
toppriority_administrativeleader <- datasetAF[,18]
controll_minister<-datasetAF[,19]
controll_administrativeleader <-datasetAF[,20]
experienced_autonomy <-datasetAF[,21]
societalvalue <- datasetAF[,22]
personalvalue <- datasetAF[,23]
chairman1<- datasetAF[,24]
projectleader1 <- datasetAF[,25]
agency1 <- datasetAF[,26]
ruimteclaim1<- datasetAF[,27]
hardeb<- datasetAF[,28]
zachteb<- datasetAF[,29]
kwaltb<- datasetAF[,30]
logistieksysteem1<- datasetAF[,31]

#attach vectors to the network netmmm#
set.vertex.attribute(netmmm,"yearsexperience",years_experience)
set.vertex.attribute(netmmm,"yearsrepresentative",years_member_AF)
set.vertex.attribute(netmmm,"controllminister",controll_minister)
set.vertex.attribute(netmmm,"controlladministrativeleader",controll_administrativeleader)
set.vertex.attribute(netmmm,"persvalue",personalvalue)
set.vertex.attribute(netmmm,"societvalue",societalvalue)
set.vertex.attribute(netmmm,"toppminister",toppriority_minister)
set.vertex.attribute(netmmm,"toppadml",toppriority_administrativeleader)
set.vertex.attribute(netmmm,"autonomy",experienced_autonomy)
set.vertex.attribute(netmmm,"chairman",chairman1) #with this vertex attribute, we identify the chairman in the network
set.vertex.attribute(netmmm,"projectleader",projectleader1)
set.vertex.attribute(netmmm,"agency",agency1) 
set.vertex.attribute(netmmm,"ruimteclaim",ruimteclaim1) 
set.vertex.attribute(netmmm,"hardebelangen",hardeb) 
set.vertex.attribute(netmmm,"zachtebelangen",zachteb) 
set.vertex.attribute(netmmm,"kwaltbelangen",kwaltb) 
set.vertex.attribute(netmmm,"logistieksysteem",logistieksysteem1) 

#add covariates beliefhomophily, trust and earlier billateral contact patterns
beliefhomophily<- as.matrix(read.table("beliefhomo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
trust<- as.matrix(read.table("trust.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
earlierbillateralcontact<- as.matrix(read.table("earliercontact.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))
actorimportance<- as.matrix(read.table("actor_impo.csv", header=TRUE, row.names=1, sep=";", stringsAsFactors=FALSE))

diag(beliefhomophily)<- 0
diag(trust)<- 0
diag(earlierbillateralcontact)<- 0
diag(actorimportance)<- 0

########################3 end of preparation dimension of constructive conflict###############################

##############################R-Analysis dimension of construction#################################

################first ERGM analysis with only the endogeneous term of 'edges'###############
c1=ergm(net~edges,
         control=control.ergm(seed=1603, force.main = T))#force-main=T wanneer geen dependence term is included.
summary(c1)
mcmc.diagnostics(c1)#
sink(file="c1.txt") 
print(summary(c1))
c1$coef
print(c1$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(c1)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("c1.pdf")
gof_c1 <- gof(c1 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_c1)
print(mcmc.diagnostics(c1))
dev.off()

setEPS()
postscript("c1.eps")
par(mfrow = c(3, 2))
plot(gof_c1)
dev.off()
dev.off()

win.metafile("c1.wmf")
par(mfrow = c(3, 2))
plot(gof_c1)
dev.off()
dev.off()

#edges significant means, according to Harris (2014:46), "..in this case trhe coeffcient for the edges term is negative and signifcant, indicating that the densitiy of the network is below 50%".

####################Second ERGM analysis with the other endogeneous terms (edges, mutual, and gwodegree)##############
c2=ergm(net~edges+gwodegree(0.7)+mutual,
        control=control.ergm(seed=1603, force.main = T))
summary(c2)
##The gwodegree(decay) term also includes the normal gwodegree in the analysis.
mcmc.diagnostics(c2)#
sink(file="c2.txt") 
print(summary(c2))
c2$coef
print(c2$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(c2)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("c2.pdf")
gof_c2 <- gof(c2 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_c2)
print(mcmc.diagnostics(c2))
dev.off()

setEPS()
postscript("c2.eps")
par(mfrow = c(3, 2))
plot(gof_c2)
dev.off()
dev.off()

win.metafile("c2.wmf")
par(mfrow = c(3, 2))
plot(gof_c2)
dev.off()
dev.off()

############################Third ERGM analysis with nodal covariate effects############################# 

c3=ergm(net~edges+mutual+gwodegree(0.7)+ttriple+nodecov("agency")+nodecov("yearsrepresentative")+nodecov("toppminister")+nodecov("persvalue")+nodecov("chairman"),
        control=control.ergm(seed=1603, force.main = T))
summary(c3)
mcmc.diagnostics(c3)#
sink(file="c3.txt") 
print(summary(c3))
c3$coef
print(c3$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(c3)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("c3.pdf")
gof_c3 <- gof(c3 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_c3)
print(mcmc.diagnostics(c3))
dev.off()

setEPS()
postscript("c3.eps")
par(mfrow = c(3, 2))
plot(gof_c3)
dev.off()
dev.off()

win.metafile("c3.wmf")
par(mfrow = c(3, 2))
plot(gof_c3)
dev.off()
dev.off()

################Fourth ERGM analysis with dyadic and nodal covariate effects##################
c4=ergm(net~edges+mutual+gwodegree(0.7)+ttripple+nodecov("chairman")+edgecov(actorimportance)+nodecov("yearsrepresentative")+edgecov(beliefhomophily)+edgecov(trust),
        control=control.ergm(seed=1603, force.main = T))

summary(c4)
mcmc.diagnostics(c4)#
sink(file="c4.txt") 
print(summary(c4))
c4$coef
print(c4$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(c4)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("c4.pdf")
gof_c4 <- gof(c4 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_c4)
print(mcmc.diagnostics(c4))
dev.off()

setEPS()
postscript("c4.eps")
par(mfrow = c(3, 2))
plot(gof_c4)
dev.off()
dev.off()

win.metafile("c4.wmf")
par(mfrow = c(3, 2))
plot(gof_c4)
dev.off()
dev.off()

#####################Fifth ERGM analysis - Optimized dyadic covariate effects model - including the term 'in-degreepopularity'################
###In comparison to model 4, there is room for improvement. This room for improvement is, especially, visible in the higher values of the in-degree parameter. Therefore, I decide to control for these higher values, and therewith optimize the model to better capture the emergent learning interaction patterns in the FSSPP network regarding the socio-cognitive behaviour of construction, by including the endogenous network statistic of in-degree popularity in the analysis.
###In-degree popularity is an endogenous network statistic which assumes that nodes with a higher in-degree are more attractive for others to send a tie to. In our analysis, this would imply that representatives that receive a lot of learning ties from network alters, are more attractive for other network alters to also send a learning tie to.
c5=ergm(net~edges+mutual+idegreepopularity+gwodegree(0.7)+ttriple+nodecov("yearsrepresentative")+nodefactor("chairman", base=1)+edgecov(beliefhomophily)+edgecov(trust),
        control=control.ergm(seed=1603, force.main = T))
summary(c5)
mcmc.diagnostics(c5)#
sink(file="c5.txt") 
print(summary(c5))
c5$coef
print(c5$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(c5)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("c5.pdf")
gof_c5 <- gof(c5 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_c5)
print(mcmc.diagnostics(c5))
dev.off()

setEPS()
postscript("c5.eps")
par(mfrow = c(3, 2))
plot(gof_c5)
dev.off()
dev.off()

win.metafile("c5.wmf")
par(mfrow = c(3, 2))
plot(gof_c5)
dev.off()
dev.off()

#####The optimized dyadic covariate model (model 5) is the model that is presented in the main text of the article as this is the model that best captures the interactive dynamics in the FSSPP network regarding the socio-cognitive behaviour of 'construction' between the network members. 


#####################################ERGM Analyses of dimension co-construction####################################################
#########################first ERGM analysis with the endogeneous terms 'edges'##########################
cc1=ergm(netmm~edges, 
          control=control.ergm(seed=1603, force.main = T))
summary(cc1)
mcmc.diagnostics(cc1)#
sink(file="cc1.txt") 
print(summary(cc1))
cc1$coef
print(cc1$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cc1)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cc1.pdf")
gof_cc1 <- gof(cc1 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cc1)
print(mcmc.diagnostics(cc1))
dev.off()

setEPS()
postscript("cc1.eps")
par(mfrow = c(3, 2))
plot(gof_cc1)
dev.off()
dev.off()

win.metafile("cc1.wmf")
par(mfrow = c(3, 2))
plot(gof_cc1)
dev.off()
dev.off()

#########################Second ERGM analysis with the other endogeneous terms (edges, mutual, and ttriple)##########################
cc2=ergm(netmm~edges+ttriple+mutual,
             control=control.ergm(seed=1603, force.main = T))
summary(cc2)
mcmc.diagnostics(cc2)#
sink(file="cc2.txt") 
print(summary(cc2))
cc2$coef
print(cc2$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cc2)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cc2.pdf")
gof_cc2 <- gof(cc2 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cc2)
print(mcmc.diagnostics(cc2))
dev.off()

setEPS()
postscript("cc2.eps")
par(mfrow = c(3, 2))
plot(gof_cc2)
dev.off()
dev.off()

win.metafile("cc2.wmf")
par(mfrow = c(3, 2))
plot(gof_cc2)
dev.off()
dev.off()

#########################Third ERGM analysis  with nodal covariate effects##########################
cc3=ergm(netmm~edges+mutual+ttriple+nodefactor("chairman", base=1)+nodecov("agency")+nodecov("persvalue")+nodecov("toppminister")+nodecov("yearsrepresentative"),
         control=control.ergm(seed=1603, force.main = T))
summary(cc3)
mcmc.diagnostics(cc3)#
sink(file="cc3.txt") 
print(summary(cc3))
cc3$coef
print(cc3$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cc3)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cc3.pdf")
gof_cc3 <- gof(cc3 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cc3)
print(mcmc.diagnostics(cc3))
dev.off()

setEPS()
postscript("cc3.eps")
par(mfrow = c(3, 2))
plot(gof_cc3)
dev.off()
dev.off()

win.metafile("cc3.wmf")
par(mfrow = c(3, 2))
plot(gof_cc3)
dev.off()
dev.off()

################Fourth ERGM analysis with dyadic and nodal covariate effects##################

cc4=ergm(netmm~edges+edgecov(beliefhomophily)+edgecov(actorimportance)+edgecov(trust)+mutual+ttriple+nodefactor("chairman", base=1)+nodecov("agency")+nodecov("persvalue")+nodecov("toppminister")+nodecov("yearsrepresentative"),
         control=control.ergm(seed=1603, force.main = T)) 
summary(cc4)
mcmc.diagnostics(cc4)
sink(file="cc4.txt") 
print(summary(cc4))
cc4$coef
print(cc4$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cc4)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cc4.pdf")
gof_cc4 <- gof(cc4 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cc4)
print(mcmc.diagnostics(cc4))
dev.off()

setEPS()
postscript("cc4.eps")
par(mfrow = c(3, 2))
plot(gof_cc4)
dev.off()
dev.off()

win.metafile("cc4.wmf")
par(mfrow = c(3, 2))
plot(gof_cc4)
dev.off()
dev.off()

#####Model 3 is the model that is presented in the main text of the article as this is the model that best captures the interactive dynamics in the FSSPP network regarding the socio-cognitive behaviour of 'co-construction' between the network members. 


###############################ERGM Analyses dimension constructive-conflict###############################
############First ERGM analysis with only edges as term#####################

cconf1=ergm(netmmm~edges, 
         control=control.ergm(seed=1603, force.main = T))
summary(cconf1)
mcmc.diagnostics(cconf1)#
sink(file="cconf1.txt") 
print(summary(cconf1))
cconf1$coef
print(cconf1$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cconf1)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cconf1.pdf")
gof_cconf1 <- gof(cconf1 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cconf1)
print(mcmc.diagnostics(cconf1))
dev.off()

setEPS()
postscript("cconf1.eps")
par(mfrow = c(3, 2))
plot(gof_cconf1)
dev.off()
dev.off()

win.metafile("cconf1.wmf")
par(mfrow = c(3, 2))
plot(gof_cconf1)
dev.off()
dev.off()

################Second analysis including other endogeneous terms##################
cconf2=ergm(netmmm~edges+mutual+gwodegree(0.7), 
         control=control.ergm(seed=1603, force.main = T))
summary(cconf2)
mcmc.diagnostics(cconf2)#
sink(file="cconf2.txt") 
print(summary(cconf2))
cconf2$coef
print(cconf2$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cconf2)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cconf2.pdf")
gof_cconf2 <- gof(cconf2 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cconf2)
print(mcmc.diagnostics(cconf2))
dev.off()

setEPS()
postscript("cconf2.eps")
par(mfrow = c(3, 2))
plot(gof_cconf2)
dev.off()
dev.off()

win.metafile("cconf2.wmf")
par(mfrow = c(3, 2))
plot(gof_cconf2)
dev.off()
dev.off()

###########Third analysis with exogenous terms##############################
cconf3=ergm(netmmm~edges+mutual+gwodegree(0.7)+nodecov("toppminister")+nodecov("persvalue")+nodecov("agency")+nodecov("chairman")+nodecov("yearsrepresentative"), 
            control=control.ergm(seed=1603, force.main = T))
summary(cconf3)
mcmc.diagnostics(cconf3)
sink(file="cconf3.txt") 
print(summary(cconf3))
cconf3$coef
print(cconf3$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cconf3)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cconf3.pdf")
gof_cconf3 <- gof(cconf3 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cconf3)
print(mcmc.diagnostics(cconf3))
dev.off()

setEPS()
postscript("cconf3.eps")
par(mfrow = c(3, 2))
plot(gof_cconf3)
dev.off()
dev.off()

win.metafile("cconf3.wmf")
par(mfrow = c(3, 2))
plot(gof_cconf3)
dev.off()
dev.off()


################Fourth ERGM analysis with dyadic and nodal covariate effects##################
cconf4=ergm(netmmm~edges+mutual+gwodegree(0.7)+nodecov("chairman")+edgecov(actorimportance)+nodecov("yearsrepresentative")+edgecov(beliefhomophily)+edgecov(trust), 
             control=control.ergm(seed=1603, force.main = T))
summary(cconf4)
mcmc.diagnostics(cconf4)#
sink(file="cconf4.txt") 
print(summary(cconf4))
cconf4$coef
print(cconf4$mle.lik) # also print the log likelihood
print(mcmc.diagnostics(cconf4)) #also print diagnostics
sink()

#goodness of fit statistics in various file formats.
pdf("cconf4.pdf")
gof_cconf4 <- gof(cconf4 ~ idegree + odegree + espartners + dspartners + distance, verbose = FALSE)
par(mfrow = c(3, 2))
plot(gof_cconf4)
print(mcmc.diagnostics(cconf4))
dev.off()

setEPS()
postscript("cconf4.eps")
par(mfrow = c(3, 2))
plot(gof_cconf4)
dev.off()
dev.off()

win.metafile("cconf4.wmf")
par(mfrow = c(3, 2))
plot(gof_cconf4)
dev.off()
dev.off()

#####Model 4 is the model that is presented in the main text of the article as this is the model that best captures the interactive dynamics in the FSSPP network regarding the socio-cognitive behaviour of 'co-construction' between the network members. 
